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Power law generalized covariance functions provide a simple model for describing the local 
behavior of an isotropic random field. This work seeks to extend this class of covariance functions 
to spatial-temporal processes for which the degree of smoothness in space and in time may 
differ while maintaining other desirable properties for the covariance functions, including the 
availability of explicit convergent and asymptotic series expansions. 
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1. Introduction 

Intrinsic random functions [1, 21] provide a popular class of models for spatial processes. 
These non-stationary random processes are specified by their generalized covariance func- 
tions, which determine the variances of certain linear combinations of the process (see 
Section 2.1 for details). A particularly simple and, therefore, useful class of generalized 
covariance functions is the power law class, for which the "covariance" between two obser- 
vations is proportional to the Euclidean distance between the points raised to some power 
(unless the power is an even integer). More specifically, for x > indicating intcrpoint 
distance, C, > and N the set of positive integers, the function 7^ 

fr(-C).T 2C , c ^n, 

l C (x) = { 2(-l)W 3Cl (1) 
y £j log x, C G N 

gives a valid generalized covariance function in any number of dimensions. Despite its 
simplicity, this class of models has the important virtue of admitting a broad range of 
local behaviors for the process, which is critical when, for example, considering proper- 
ties of spatial interpolants [30]. Specifically, the larger the value of £, the smoother the 
process, so that, for example, the process is m times mean square differentiable in any 
direction if and only if £ > m. A generalized covariance function can be written as the 
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Fourier transform of a positive measure; although, in contrast to the stationary setting, 
the measure might not have finite total mass. In particular, in d dimensions, the mea- 
sure corresponding to 7^ has density with respect to Lebesgue measure proportional to 
|u>| _2< ' , where u> £ M. d is the spatial frequency. 

The goal of this paper is to find a good extension of 7^ to the space-time setting. The 
problem is made difficult by what I will mean by "good." The first requirement is that 
the class of models includes members allowing any degree of smoothness in space and any 
(possibly different) degree of smoothness in time. Specifically, for G(x, t) a generalized 
covariance function of spatial lag x and temporal lag t, and any positive Ci)C2>Cx and 
C2, the class of generalized covariance functions should include a member that satisfies 

G(x,0) = Ci 7Cl (|x|) (2) 

and 

G{Q 1 t) = C 2lC2 {\t\) (3) 

or, failing that, that (2) holds asymptotically as x — > 0, and (3) holds asymptotically as 
*->0. 

The second requirement is that G be smoother away from the origin than it is 
at the origin. This requirement is needed to avoid the kinds of anomalies described 
in [31] for covariance functions that are not smoother away from the origin. As a 
simple example of this kind of anomaly, consider the covariance function onlxl 
given by K(x,t) = exp(— \x\ — \t\). Write corr for correlation, and define p(x,t) = 
lim e ^o corr{Z(0, e) — Z(0, 0), Z(x, t + e) — Z(x,t)}. Then for x 7^ 0, straightforward cal- 
culations show p(x, t) = for t 7^ and p(x, 0) = c - ' 1 ' . The discontinuity in this limiting- 
correlation is due to the fact that K(x,t) has a similar discontinuity in its first derivative 
in the t direction everywhere along the x axis as it does at the origin. I consider such 
a discontinuity in p as unrealistic for most natural space-time processes. In particular, 
Stein and Handcock [32] and Stein [29] give examples showing how this lack of continuity 
away from the origin can lead to optimal (kriging) predictors with undesirable properties. 
All separable space-time covariance functions, that is, those that factor into a function of 
space and a function of time such as c - '^ - '*', have a similar problem unless the process 
is infinitely diffcrentiable [31], page 311. Furthermore, many non-separable space-time 
covariance functions proposed in the literature share this problem [31], pages 311-312. 
For a space-time covariance function with different degrees of smoothness in space and 
time, that is, satisfying (2) and (3) with £1 7^ £2, it is not so clear what one should mean 
by the function being smoother away from the origin than at the origin. This issue is 
addressed in Section 2.4. 

The smoothness of a covariance function away from the origin is closely related to 
regularity properties of the corresponding spectral density at high frequencies [31]. In- 
deed, Stein [29] argues that the spectral domain provides a more natural approach for 
considering the appropriateness of various models for space-time covariance functions. 
Specifically, [29] gives the following condition as a plausible requirement for the spectral 
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density f{u)) of a natural process in space or space-time: for every R < oo, 

/(w + 1/) 



lim sup 



/(«) 



- 1 



= 0. (4) 



That is, / changes slowly (on a relative scale) at high frequencies. This condition excludes, 
for example, separable space-time models. 

The third and final requirement for the generalized covariance functions is that they 
can be computed accurately and efficiently via, for example, series expansions, to allow 
them to be applied routinely to large space-time datasets. In particular, representations 
of the function as an integral will not be considered an adequate solution to the problem. 
As best as I am aware, no existing class of generalized covariance functions satisfies all 
three of these requirements. 

It will be convenient to avoid explicitly distinguishing between space and time and 
consider processes Z(x,y) on R dl x R d2 for positive integers d\ and cfe- Stein [31] proposed 
the following class of spectral densities (with a different parameterization) as a flexible 
parametric model for stationary space-time processes: 

for <7i,<72 positive, ftl + /3| > and 2v > d\/a\ + (fe/a^; this last condition being nec- 
essary and sufficient (given the positivity constraints on the other parameters) for q to 
be integrable. The parameters fix and /?2 are inverse range parameters, o\ and o~i are 
scale parameters and ct\ , ai and v together describe the smoothness of the process in 
x and y. More specifically, the process is p times mean square diffcrentiablc in each 
component of x if and only if 2v > (d\ + 2p)/a\ + dija.2 and, similarly, is p times mean 
square diffcrentiablc in components of y if and only if 2v > d\ju\ + (^2 + 2p)/c*2 [31]. 
Furthermore, when a± and ui arc integers, the resulting covariance function is infinitely 
diffcrentiablc away from the origin ([31], Proposition 4, although this result also follows 
from [28], Theorem 1.1). All models in the class (5) satisfy the spectral condition (4). 
Porcu [22] describes more general approaches to obtaining valid spectral densities that 
could be useful in the space-time context. 

An obstacle to using the class of models (5) is the lack of an explicit expression for 
the corresponding covariance functions. Except when a\ = Q'2 = 1 and some special cases 
with v an integer, a.\ = 1 and ai = 2, I am unaware of any cases for which an explicit 
expression has been written down [31]. For rational spectral densities, which includes the 
model (5) when ai,«2 and v are all integers as a special case, the covariance function 
can be expressed as the solution of a certain set of equations that reduce to a partial 
differential equation when the rational function is just a reciprocal of a polynomial [25, 
26]. However, writing down an explicit solution for all integers a\,a2 and v is not a 
simple task. Ma [18] gives explicit expressions for the covariance functions of space- 
time processes with rational spectral densities in some limited special cases that do 
not include any instances of (5). For all of the cases treated in [18], the covariance 



4 



M.L. Stein 



functions are not smoother away from the origin than at the origin in the sense defined 
in Section 2.3, and the spectral densities do not satisfy (4). Even if one had explicit 
expressions for the covariance functions of all rational spectral densities, these covariance 
functions (asymptotically in a neighborhood of the origin) satisfy (2) and (3) only for a 
countable nowhere dense set of (Ci : C2) values. Kclbcrt, Lconcnko and Ruiz-Medina [12] 
consider the case d± = d>2 = 1 ; ot\ = 1, 0:2 = 2 and \}\ = as a stochastic fractional heat 
equation in some detail, but they only give integral representations for the covariance 
functions. They consider additionally setting /?2 = and note that the resulting random 
field has a self-similarity property when v € (f,|)- Christakos [3], page 225, mentions 
the case with d\ = l,ai = v = 1 and j3\ = as a possible model for spatial-temporal 
processes, but derives no results for these models. 

The parameters /?i and /?2 are range parameters that do not affect the local behavior 
of the process, so consider setting /?i = P2 — in (5), yielding 



These spectral densities are not integrable in a neighborhood of the origin, so they corre- 
spond to generalized covariance functions. Section 2.1 gives some background on gener- 
alized covariance functions. Theorem 1 in Section 2.2 gives a convergent power series for 
the generalized covariance function corresponding to (6) when a.% = 1, a± > 1 and y 7^ 
and separate explicit formulae to cover the case y = 0. For completeness, the known 
result [6] for the generalized covariance function when a\ = 1 is also given in Theorem 1. 

Section 2.3 shows how the generalized covariance functions in Theorem 1 can, in most 
cases, be written in terms of the H- function [14], a generalization of the generalized 
hypergeometric function that is sometimes called Fox's i/-function [5]. This result is used 
to obtain asymptotic series for these generalized covariance functions and to motivate 
a conjecture extending Theorem 1 to the case a\ < 1. Section 2.4 shows that for any 
positive Ci,C2,Ci and C2, one can find a spectral density of this form satisfying (2) 
and (3). Furthermore, in a sense made precise in Section 2.4, the resulting generalized 
covariance function is shown to be smoother away from the origin than at the origin. 
Finally, these results are used to show that at least some of these covariance functions 
avoid what [13] calls the "dimple" that occurs in some proposed space-time covariance 
functions, which is a lack of monotonicity in the correlation structure that might often 
be viewed as unnatural. 

Section 3 discusses some limitations and possible extensions of the generalized covari- 
ance functions considered in Section 2. This section also touches on some of the difficulties 
in using the series expansions to compute these functions quickly and accurately. Sec- 
tion A.l gives a proof of Theorem 1, and Section A. 2 collects some needed material on 
ii/- functions. 

There is a substantial recent literature on the development of space-time covariance 
functions with explicit representations in terms of well-known special functions. Some 
references include [3, 4, 7, 9, 10, 15-19, 23, 24, 27, 31, 33]. Although these works consider 
a broad range of models for space-time covariance functions, none of them give a class 
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of covariance functions meeting the criteria set forth in this section of (locally) satisfying 
(2) and (3) for all positive £i and £2 as well as being smoother away from the origin 
than at the origin. As noted in [31], perhaps [16] comes closest to this goal, in that this 
paper gives a class of models with d% = \ including elements satisfying, asymptotically in 
a neighborhood of the origin, (2) with Ci = I an d (3) with < C2 < \ , and the covariance 
functions are infinitely diffcrentiablc away from the origin. 



2. Theoretical results 

2.1. Generalized space— time covariance functions 

Intrinsic random functions and generalized covariance functions have been a standard 
tool in geostatistics since Matheron's pioneering work [21]. These processes are nearly 
stationary in the sense that variances of some class of linear combinations of the process 
are translationally invariant. Specifically, for a random field Z on R d (so that d = d\ + di 
is the total number of dimensions for space-time processes), call XjZ(zj) an autho- 

rized linear combination of order fc, or ALC-fc, if XjP(zj) = for every polynomial 

P of order at most fc. A function G on M. d is called a generalized covariance function 
of order fc, or GC-fc, if, for every ALC-fc, X)™j=i XeXjG(z^ — zj) > 0. A process Z for 
which Var{£™ =1 XjZ(zj)} = E" J=1 X e XjG(z e - Zj) > and £{£" =1 X j Z{z j )} = for 
every ALC-fc is said to be an intrinsic random function of order fc, or IRF-fc, with G 
as its GC-fc. GC-fcs are not unique; if G is a GC-fc for Z, then so is G plus any even 
polynomial of degree 2fc in x. A function / is the spectral density of a GC-fc G if 



n ,. 

V XgXjGize-Zj) = I 



n 

3=1 



/(w)dw 



for every ALC-fc ^2" = iXjZ(zj). A nonnegative even function / is the spectral density 
for a real-valued GC-fc if and only if /(w)|w| 2fe + 2 /(l + |w| 2fe + 2 ) is integrable [21]. Of 
course, if Z is an IRF-fc, it is also an IRF-fc' for all integers fc' > fc. Let [x\ indicate 
the greatest integer less than or equal to x. For the spectral density (6) and i = 1,2, 
define fej = \a.i{v — d\ / (2a\) — d%j '(2a,2)}\ and fco = max(fci, k%). Straightforward calculus 
shows that / in (6) satisfies the integrability condition for an IRF-fc if and only if fc > fco- 
Christakos [2] considers an extension of the notion of IRFs to the space-time setting in 
which one essentially allows a different degree of differencing in space and in time, but 
this concept is not used here. 



2.2. Main theorem 

This section gives a series expansion for the generalized covariance function corresponding 
to (6) when U2 = 1, ol\ > 1 and ai = 02 = 1- Extending the result to other positive values 
of a 1 and 02 is trivial. Define the function M. v (y) = y u JC^(y), where K, v is a modified 
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Bessel function. When v > 0, M„ is often called the Matern covariance function with 
smoothness parameter v. Set 9 = v — ^d,2 and 9' = 9 — di/(2ai) so that kg = \_ai6'\ is 
the (minimal) order of the IRF corresponding to this spectral density. As noted in [31] 
for the more general model (5), when a 2 = 1, the Fourier transform with respect to u> 
can be carried out explicitly. Specifically, for t ^ 0, 

|r| 1 + \uj\) e y da;= 2 e- lT ^ T \ 2ai e ■ ( 7 > 
Define rtj = |jq — Xj |, S£j = \yt — y 3 - 1 and, for t > 0, 

(-(l/4)t 2 ) m 



A d ( i ) = 2 ( d - 2 )/ 2 r(iA-( rf - 2 )/ 2 j (rf _ 2)/2 (t) = f: 

\ / ™ _ f 



- 77l!(d/2) r 
m— 



(8) 



where, for any real a and positive integer j, (a)j = a(a + 1) • • • (a + j — 1) and (a)o = 1. 
For every ALC-fco Ylj=i ^jZ{x-j,yj), the GC-fco G corresponding to the spectral density 
(|t| 2qi + \u>\ 2 )-" satisfies 

n 

AMj G(x £ - Xj , y e - y 3 ) 



f f V AfAj-c^'^-^'+^'^-^'drl 2 " 1 + |w| 2 )- 1/ du;dx 



2»" 1 r(i/) 7 Rdl ^ J |r| 2 ^ 

POO 71 

Jo « •_, 



(9) 

47t (d 1 +d 2 )/2 



2*r(*/)r((i/2)d 1 ) 7 

47T (d 1 +d 2 )/2 
= 2 e r( i /)r((l/2)d 1 )a 1 

x / V \i\jt dl l ai ~ 2e ~ 1 kd l (rejt 1 / ai )Me(stjt) dt, 

where the second step uses (7), the third basic results on Fourier transforms of isotropic 
functions [30], Section 2.10, and the last step the change of variables t = u ai . What one 
would like to do is substitute (8) into (9) and integrate termwise, but justifying this 
interchange requires considerable care. 
Define 

. _ 7t( dl+d2 )/ 2 r((d 1 + 2m)/(2ai)) 
Cm(ai) ~ m!r(m + (l/2)di) 
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and the digamma function tp by ip(z) = ^Mogr(z). Equations (10)-(12) are proven in 
Section A.l. Equations (13) and (14) are taken from [6], Chapter II, Section 3.3, equations 
(2) and (11). 

Theorem 1. For the spectral density /(t,u>) = (|x| 2qi + \u)\ 2 )~ u on R dl x R d2 with 
on > 1, 6' = v — d\j (2ai) — c?2/2 > and kg = \_a\9' J , a corresponding GC-ka is given by 
G(\x\, |y|) for the function G, defined by the following equations. First consider otx > 1. 
For s > 0, 

Cm(ai) J / 1 \ /I 



m=0 

with 7(-) defined by (1). When a\9' > fco 



"5m>=E^H?J j — (?) (10) 



n (d 1+ d 2 )/2 T(g \ 

G (r, o) = w zr- r«io> l^r), (11) 



r(i/)r(aifl) 

and when a\9' = ko, 



r(i/ ) I V 2 / J I V 2 / «i V 2 "i 



Finally, when ol\ = 1 and 0' is noi an integer, 



+ — ^(1) - V ( fco + ^di ) - V^fco + 1) 



7T 



i+d 2 +2)/2 



and when ot\ — 1 and 0' is an integer, 



(12) 



^ S ) = - sin W V + 1)^)2^ ^ + 82)6 > (13) 



C_-| \fc+l„(d 1 +d2)/2 

G(r,s)= { - l) e ,^ v)22S , (r* +S y log(r 2 +.s 2 ). (14) 



2.3. iJ-functions 

The function G can, in most cases, be written in terms of H- functions [14, 20]. Sec- 
tion A. 2 gives the definition of //-functions as a contour integral and some other needed 
information about the functions. If 

m 

t-\ 7^ 9' for all whole numbers £ and m, (15) 
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then (45) in Section A. 2 is satisfied, and 

/ 



H. 



2.2 



V 



(1,1), 1^,1 



di 1 

2«i ' ai / 



1 



is well defined. Since, by (47), A = — 2, which is negative for a± > 1, [14], Theorem 1.4, 



applies, yielding, for z ^ 0, 
( 



H. 



2.2 



V 



(1,1), ('$-,1 



dl 1 

2ai ' ai / ' 



where 



/j lfc = (-l) fe air 



di +2fc 
2ai 



1 

ai 



r( -#' + — 



(16) 



Comparing this result to (10) yields (assuming (15)) 

/ 



G(r,s) 



n (d 1+ d 2 )/2 ( X ^28' 

I»ai 



r 



H 



2.2 



V 



((l/2)s) 2/ ° 
((l/2)r) 2 



(1,1), ('£,1 



2ai ct\ ) ' 



1 



(17) 



/ 



where this result holds for r = by continuity. For the parameter values of the H- 
function in (17), from (46) and (47) in Section 3.2, A < and a* > 0, so that by [14], 
Theorem 1.11, Hl^\{z) has an asymptotic expansion as z — > for | argz| < ^a*n. To avoid 
complications, let No be the set of nonnegative integers, assume 8 ^ No and that for all 
I G No, (6' — t)a x i N and (0 - e)a t £ N. Then, as z -> 0, 



H. 



2.1 



/ 



V 



(1,1), 



,1 



( d t 




-"'-) 


V2"i 




a x J 



Y,{h* u z lai+dl/2 + K e z^- e '^}, (18) 



where, using the duplication formula for T, 

{-lYa x Y{-e-i)T{{l/2)d 1 +la l ) 



h* 



tW{-tai) 

ai sin(7t£ai)r((l/2)d 1 + l ai )T{lai + 1) 
sm(nO)T{9 + l + l)£\ 



(19) 
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and 

_ (-iy ai r(8-i)r((£-6>)a 1 ) 

2e eT(ai(9-t)) 

_ on sin{7Tai(6» - £)}F((t - 0')a x )V{{l - 0)on + 1) 
~~ sin(7t6»)r(l -e + ey\ ' 

From (19), /i* = and, if ol\ is an integer, h* u = for all I. For 9 G No, [14], (1.8.2) applies, 
yielding an asymptotic expansion with logarithmic terms, but the result is rather messy 
and is omitted here. 

When (I s ) 2 /«i/(I r ) 2 is small, (10) converges slowly and is numerically unstable. 
Specifically, for any fixed m > ot\9' and r > 0, as s J, 0, the mth term in (10) tends 
to ±oo, even though lim^o G(r, s) —> G{r, 0), which is finite. Thus, there must be a near 
canceling of large terms of opposite signs in (10) for s small, so that high precision arith- 
metic would be needed to obtain accurate results for s sufficiently small. Fortunately, the 
asymptotic expansion (18) can be used to approximate G(r,s) for s small. Substituting 
(18) into (17) and considering r > fixed, 

as s |. 0. Since h\ = 0, when 9 is not an integer, 

L0J+1 

G(r,s) = A lS 2l + B le+ i{s)+o{ S w ) (21) 
e=o 

as s 4- for some constants (depending on r ) , . . . , A + 1 , B , where B = if and only 
if Qi is an integer. The asymptotic expansion (21) also holds when 9 is an integer by [14], 
(1.8.2) and (1.4.5). 

It is apparent that (10) does not give a valid power series expansion for G when a.\ < 1, 
since it is easy to show that the individual terms in the sum do not tend to as m — » oo 
for fixed and positive r and s. Nevertheless, the representation in terms of H- functions 
given by (17) may still be valid for a± < 1 when 2a\9 > d\ (so that 0' > 0). Excluding 
values of (a\,v) for which (15) is not satisfied, a natural conjecture is that (17) holds 
when eti < 1 if 2a\9 > d\. A plausible approach to proving this conjecture would be to 
use analytic continuation, but this would require at the least extending the definition of 
iJ-functions to a strip of complex values of a\ containing the positive real axis. Note that 
when oil < 1, (47) implies A > and, by [14], Theorem 1.3, (18) becomes a convergent 
power series for H^- On the other hand, now (16) is no longer a convergent power 
series, but, with a* defined as in (46), it is a valid asymptotic expansion as z — > oo for 
| argz| < ±a*7t by [14], Theorem 1.7. 
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For oti = 1, (17) can be directly verified when (15) holds. Specifically, for ol\ = 1, 
using [14], Property 2.2, the right-hand side of (17) reduces to 



n (d 1 +d 2 )/2 / x 



I» 



26 ' is (s 2 (1,1) 

^ { r-2 r-0',1) 



(22) 



The parameter A defined in (47) equals and, thus, one can show that the f/-function 
has a convergent power series given by [14], Theorem 1.3, when s <r and by [14], Theo- 
rem 1.4, when s > r. Consider s <r and 9' ^ Nq. Then straightforward calculations yield 
that (22) equals 



7T 



72 



I»2 2e ' 



(-l) e T(-0' + £) (s 



2C 



71 



(di+d 2 +2)/2 



sin(7t6»')r(6l' + l)r(i/)2 2e 



which, using the binomial series, equals (13). A similar argument shows (22) equals (13) 
for s > r, and it additionally holds for s = r by continuity. 



2.4. Consequences 

One goal of this paper was to find a class of generalized covariancc functions that has a 
member satisfying (2) and (3) for all Gi,G2,Ci and C2 positive. In fact, the functions of 
the form G(&ix, &2y) with u\ > 1 and b\ and 62 positive satisfy this requirement. To prove 
this, first suppose (1 ^ (2 and, without loss of generality, take Ci > (2- Theorem 1 gives 
explicit expressions for positive constants D\ and D2 such that G(x, 0) = -Di7 ai e'(i|x|) 
and G(0,y) = -D27e'(2"lyD- Setting 8' = (2 and ot\ = C1/C2 achieves the desired degree of 
smoothness in all directions. From (10), G(0,&2s) = D2(^b2)'' 2 ^( 2 (s), so by appropriate 
choice of &2, one obtains G(0,&2s) = G27<; 2 (s). For £1 = u\& not an integer, by (11), 
G(6ir, 0) = Di^bi^j^ (r), so by appropriate choice of bi, one gets G(&ir, 0) = Gi7^ (r). 
When (1 = a\Q' is an integer, by (12), G{b\r, 0) = Z?i(|6i) 2 ^ 1 7^ 1 (r) plus some constant 
times r 2 '' 1 , so that b\ can be chosen to make G(6ir, 0) = Gi7f 1 (r) plus some constant 
times r 2 ^ 1 . When Ci = C2, set ol\ = a.2 = 1 and £i = v — \d\ — By (13), there exists 
D > such that G(&iX, b 2 y) = Dj Cl (^|6ix| 2 + |6 2 y| 2 ). As before, one can clearly choose 
b\ and 62 so that (2) and (3) are satisfied as long as one ignores an even polynomial of 
degree 2£i when d is an integer. Since a GC-fc is only identified up to even polynomials 
of degree at most 2k, it is fair to say that the class of generalized covariance functions 
corresponding to (6) with 012 = 1 and ot\ > 1 includes members satisfying (2) and (3) for 
all Ci >C 2 ,Gi and G 2 . 

Now consider in what sense members of G achieve the goal, identified in Section 1, 
of being smoother away from the origin than they are at the origin. If a function is 
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not infinitely differentiable in any direction at the origin but is infinitely differentiable 
everywhere but the origin, then one might say without controversy that such a function is 
smoother away from the origin than at the origin. Thus, when a\ £ N, the issue is settled. 
But when ot\ is not an integer, B in (21) is not 0, and G is not infinitely differentiable 
away from the origin. 

One way to describe the smoothness of a function is by its pointwisc Holder exponent 
s. Consider a function / from R d to R, s > and x e R d . Then / £ C s (x ) if and only if 
there exists e > and a polynomial P of degree less than [s\ and a constant C such that 
|/(x) — P(x — x )| < C|x — x | s for all x satisfying |x — x | < e. The Holder exponent 
of / at xq, which I will denote by HE(xo,/), equals sup{s: s £ C s (xo)}. Because the 
degree of smoothness at the origin of G varies in different directions, it will not suffice 
to compare the Holder exponent at the origin to the Holder exponent elsewhere. To 
avoid the kind of anomaly described in Section 1, consider the smoothness of G in each 
direction separately. Specifically, for vectors Zo,Zi <E M. dl+d ' 2 and (generalized) covariance 
function K, consider the function of t £ R given by K(t; Zq, zi) = K(zq + tz\). Define 
to be a vector of zeroes whose length is apparent from context. Then I claim that, in the 
present setting, a useful notion of K being smoother away from the origin than at the 
origin is 

HE(0,if(-;z o ,Zi)) >HE(0,K(-;0,zi)) for all z + 0,z x ^ 0. (23) 

To see why this definition might be appropriate here, consider the following gener- 
alization of the example in the introduction. Suppose K is a covariance function or a 
generalized covariance function of order and, for nonzero Zo and Z\ and < a < 2, 
X(tzi) = C + Ci|t| Q + 0(1*1") and K{z + tz x ) = D + £>i|i| a + D 2 t + o{\t\ a ) as t -> 
for some constants Co, Gi, Do, -Di and D 2 (possibly depending on zo and zi) with C\ 
and D\ nonzero. It follows that HE(0, K(-; 0, zi)) = HE(0, K(-; zo, Zi)) = a. Furthermore, 
limt_>o corr{Z(tzi) — Z(0),Z(zo + tz\) — Z(zq)} = D\/C\ ^ 0. Now suppose the lack 
of smoothness of K at Zo is localized in the sense that there exists e > such that 
HE(0, K(-;zq + (5zi,Zi)) > a for all < \6\ < e, which, as far as I am aware, holds for 
any space-time covariance function that has been proposed in the literature. This con- 
dition implies lim 4 ^ corr{Z(tzi) — Z(0),Z(z + Sz-y + tz\) — Z(z + Szi)} = for all 
< \S\ < e. Thus, when K is not smoother in the Zi direction at zo than it is at 0, 
there is a "discontinuity" in correlations of increments. If, instead, D\ = 0, which will be 
the case under (23), then this limiting correlation is for all 8 in a neighborhood of 
including 5 = 0, and no discontinuity occurs. 

Define Zj = (x.j,yj) for j = 0, 1 with Xj £ M. dl and yj £ R d2 . For G as given in Theo- 
rem 1, let <5(i;z ,zi) = G(z +tz 1 ). For z x ^ 0, by (10), HE(0,G(-;0, (xi.yi)) = 29' if 
xi ^ and HE(0, G(-; 0, (0, yi)) = 2ai0'. Now consider z ^ 0. If y ^ 0, then it is pos- 
sible to show (10) can be differentiated termwise and HE(0, G(-; (x , y ), (xi, yi))) = 
oo. Next, HE(0,G(-;(x ,0),(0,yi))) = oo by (11) or (12). Finally, if Xi ^ 0, then 
HE(0,G(-;(x ,0),(xi,yi))) > 29 + 2 by (21). Because 9 > 6', in all cases HE(0,G(-;z ,Zi)) > 
HE(0, G(-; 0, Zi)) +2 for all nonzero Zo and Zi, so (23) is more than satisfied. 

It is not clear that satisfying (23), or even the stronger condition met by G here, will 
exclude all possible "discontinuities" or other anomalies in the covariance structure, but 
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it does avoid at least the type considered here. It might be preferable to find covariance 
functions that satisfy (2) and (3) and are infinitely differentiable away from the origin, but 
for £1,(2 and (1/C2 ah irrational, I am unaware of any generalized covariance functions 
that satisfy all of these conditions. 

Next, consider the problem of the dimple in fully symmetric (even in both its argu- 
ments) stationary space-time covariance functions described in [13]. A formal definition 
of the dimple is given in [13], but the essential point is that the space-time covariance 
function K(x,t) has a dimple in, say, the time lag t if, for some fixed spatial lag x, 
K(x,t) has a local minimum in t at t = 0. This dimple implies that at the spatial lag x, 
correlation is stronger with both the near future and the near past than with the present, 
and [13] argues that such a lack of monotonicity in the covariance structure will often be 
undesirable. A dimple in the spatial lag can be defined similarly. 

For a GC-fc with k > 0, it is not clear what one should mean by a dimple, but for k = 0, 
the GC-0 that equals at the origin is just minus the semivariogram for the process: 
I Var{Z(x, t) — Z(0, 0)} = — G(x, t). If the variogram is expected to increase as one moves 
"farther away" in space-time, then GC-Os with dimples should be avoided. 

Assume fco = or, equivalcntly, a.\9' < 1, so that G is a GC-0. First, if ct\ = 1, G(r, s) is 
a decreasing function of \/r 2 + s 2 , and there is no dimple, so assume ai > 1. In addition, 
assume is not an integer so that (20) holds. To show that G(r, s) does not have a dimple, 
it suffices to show that for every s > 0, G(0,s) > G(r, s) for all r sufficiently small and, 
for every r > 0, G(r, 0) > G(r, s) for all s sufficiently small. That G(0,0) = is greater 
than G(r, 0) and G(0,s) for all positive r and s is immediate from (10) and (11). For 
fixed s > 0, from (10), 

S) " ^ S) = ^f^-^ (¥) + 0ir4) 

as r I 0. It follows from Ci(ai) > and 1/ai - 9' > that G(0,s) > G(r,s) for all r 
sufficiently small. From (20), for fixed r > 0, 

-. n . ~ . n^+ d ^/ 2 T(9-l)T((l-9') ai ) fl Y^ 6 '-^ fl \ 2 . 2 . 
G(r,0) -G(r,s) = f^^i)) UV W + ° (s ) 

as s 10. For 9 > 1, G(r, 0) — G(r, s) is clearly positive for all s sufficiently small, so 
now consider < 6 < 1. In this case, T(6 — 1) < 0, so G(r, 0) > G(r, s) follows for all s 
sufficiently small if T(a\(9 — 1)) < 0, which holds if a.\[9— 1) £ (—2m— 1,— 2m) for some 
nonncgativc integer m. This condition does not hold for all oi\ and 9' for which a\ > 1 
and a i#' < 1, but it does always hold when 9' > jqrjj-- To prove this, it suffices to show 
ai{6 - 1) > -1, which holds if & > 1 - (1 + |di)/ai- The curves 0' = 1 - (1 + |di)/ai 
and 0' = 1/ai intersect at (ai,0') = (2 + ^2, 2/(4 + di)), from which it follows that 
a±(9 — 1) > — 1 holds for all 0' < 1/ai whenever 0' > 4 _p rf . This lower bound is | for 

di = 1 and is smaller for larger di. The lower bound of | may not be too restrictive in 
practice: Brownian motion has generalized covariance function proportional to 71/2, and 
processes less smooth than Brownian motion are somewhat uncommon in applications. 
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3. Discussion 

For G as defined by (9), the series expansions (10) and (20) should, in principle, allow 
fast and accurate calculation of G, but there do not appear to be any publicly available 
programs for computing .//-functions and writing general purpose code to carry out these 
calculations would require a major effort. In particular, preliminary investigations suggest 
that considerable care needs to be taken to piece together the convergent power series 
(10) and the asymptotic expansion (20) to obtain accurate approximations for all values 
of the argument of the function. Further work would also be needed to handle those 
values for (ai,u) for which one of the expansions has a singularity or near singularity, 
including values of a.\ near 1. 

Restricting a.i = 1 in (6) was essential to the derivation of series expansions for 
the resulting generalized covariancc functions. Model (6) was, in turn, a simplifica- 
tion of (5), which includes two range parameters. Perhaps //-functions can be used 
to express, in at least some cases, the (generalized) covariance functions correspond- 
ing to these more general models. However, even the richer class of covariance func- 
tions given by (5) is inadequate for modeling many natural processes. In particular, 
these covariance functions all satisfy G(x, y) = G(x, — y) and, hence, are all what [7] 
calls fully symmetric. Any process with a predominant direction of flow will not be 
fully symmetric, so this constraint is often inappropriate. Covariance functions that 
are not fully symmetric can be generated from covariance functions that are [11, 31], 
and these approaches can, in principle, be applied to the models considered here. 
An easy extension of this model is to allow for geometric anisotropies in either x 
or y by considering G(|Ax|, |£>y|) for any d\ x d\ matrix A and any di x di ma- 
trix B. 

As noted in the Introduction, there has been quite a lot of research in recent years 
developing new classes of space-time covariance functions. In many of these works, the 
focus has been on obtaining simple closed form expressions for space-time covariance 
functions. Having closed form expressions is certainly valuable in applying the models, 
but it is critical that any such model provides a good description of the spatial-temporal 
variations of the process to which it is to be applied. Comparing various models in a 
broad range of applications is one important way to learn about which models will be 
of most use in practice, but it is also important to consider the theoretical properties 
of these models, such as their smoothness properties, both at the origin and away from 
the origin, and the presence of dimples or other possible anomalies. Finding covariance 
function models for space-time processes that allow for a different degree of smoothness 
in space and in time, possess certain other desirable properties such as (4) and are 
accurately computable using series expansions is a major challenge. The results obtained 
here perhaps provide a first step to show that it may not be necessary to sacrifice desired 
theoretical properties of space-time models in order to gain computational tractability, 
although admittedly quite a bit of work on numerical methods would be needed before 
the generalized covariance functions proposed here could be used routinely (or even not 
so routinely) in practice. 
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A.l. Proof of Theorem 1 

To prove (10) for the process Z with spectral density (6), consider the process Z\(x\, y) = 
£((a?i,0,...,0),y), which has GC-fc G 1 {x 1 ,y) = G((x u 0, . . . , 0),y) = G(\ Xl \, |y|). As- 
sume for now that fco = L a i^'J does not equal ct\Q' so that a,\& — 1 < fco < otiO' . Then 
the process Z\ is fco times mean square differcntiable in its first coordinate direction and, 
for m < fco , denote its mth mean square derivative process by Z™ (x\ , y) . The generalized 
covariance function for Z™(xi,y), denoted by G™(xi,y), can be chosen to satisfy 

GT(x 1 ,y) = (-l) m ^G 1 (x 1 ,y). (24) 

Now ^j m + +1 Gi(0,y) = for m < fco, so if one knew Gi(0, y) for m < fco, then 

Gi(xi,y) could be recovered from G 1 °(xi,y) by integration. Then, since G(x, y) has 
a version that only depends on x through |x|, one can obtain G. 

Suppose di > 1. The case d\ = 1 requires a slightly different but easier argument. 
The process Z 1 (x±,y) is an IRF-0, so its GC-0 can be taken to equal negative the 
scmivariogram of the process. Denoting |y| by s, 

-G k °(x 1 ,y)= [ f (l-e iTlXl+iu} ' y )T^ ko (\T\ 2ai + \uj\ 2 )- u dujdT 

Switching to hyperspherical coordinates with |r| =u, t\ = ucos^i and integrating over 
the angles 4>2, ■ ■ . , 4>di-i yields 

n (d 1 +d 2 -l)/2 

~ G i° (xi ' y) = 2»- 2 r( I /)r((d 1 - l)/2) 

u 2k -2 ai e'~l f cos 2fco0 lSm d 1 -2^ i ( 25) 

Jo 

x {Me(0) — cos(xiucos(/)i)Me(u ai s)} dcpidu. 

Making the change of variables a = cos</>i and using [8], 3.251.1 and 3.771.4, and the 
series expansion for the generalized hypergeometric function 1F2 yields 



cos 2k ° </>i sin" 1 - 2 0i{X e (O) - cos(x!UCOS(j)i)Me(u ai s)} d<j) 



1 



a 2k °(l - a 2 )^ dl ~ 3 ^ 2 {M g (0) - cos( Xl ua)M g (u ai s)} da 
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(26) 



1 4-1 
B[k 0+ -,^ 



{Me{0)-M e (u ai s)} 
(k + 1/2), 



E 



{k +(l/2)dM^) 



i-x\u 2 YM e {u ai s) 



where B is the beta function. 

The following properties of Me are used in the proof. For any 9 > 0, as 1 1 0, 



A4 e (i) = U r t 2r + V le (t) + o(t 2t> ) 



(27) 



r=0 



for appropriate values of the U r s and V [30], Section 2.7. In addition, for any 9 > 0, there 
exist positive constants C and I? (depending on 9) such that 



0< M e (t) <Ce~ Dt 
for all £ > 0. Furthermore, for all real 9 and all t > 0, 

d 



dt 



M fl (f) = -iMfl_i(t). 



(28) 



(29) 



Assume s > for now. The function u 2k °~ 2ai8 ' {Me(0) — Me(u ai s)} is integrable in 
u over (0, oo) since, from (28) and ko < ot\6' , it is integrable over (0, 1], and, from (27) 
and fco > ot\9' — 1, it is integrable over (l,oo). Furthermore, 



E 



(ko + 1/2)/ 



-(x 2 u 2 ) £ < cosh(xiu) 



^(*0 + (l/2)di)/(20! 
and cosh(xiu)A^e(w Ql s) is integrable over (0,oo) for ai > 1, so that 

(k + 1/2), 



E 



(-x? U 2 ) € 7W e (u Ql s) 



(A„ + (l/2)diM2*)l 



can be integrated termwise over (0, oo). Making the change of variables i = u ai , integrat- 
ing by parts and using the definition of Aig gives 



u 



2k —2a 1 8'-l 



{M e (0)~M e (u ai s)}du 
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i r°° 

= — t 2k °l^- 6 - l {Me{Q)-M e (ts)}At 
"i Jo 

n i* no 

t 2k °/ ai - 2e ' +1 M e -i(ts)dt 

i 

t (2*o-HiO/ai-»X^_ 1 ( ts )dt. 



2ai6»' - 2fc 



2ai0' - 2fc Jo 

Since (2/c + d\)/a\ — 9 > 11 — 6*1, [8], 6.561.16, applies, and one obtains 



u^-^'^iMeiO) - M e (u ai s)}du 



(30) 



U\ \OL\ ) \ 1(X\ 

Next, for i > 1, again using the change of variables t = w Ql and [8], 6.561.16, gives 
u 2ko+2i - 2aie '- l M e (u^s)du 



^0-2 / ko+£ x / 2 ko + 2e + di\ (I \ 2e '-( 2fc o+2^)/-i 



(31) 



r '^L^-e' r 



For y^O, using (25), (26), (30) and (31), 
G\°{x x ,y) 

Tx di + d i- 1 r((k + i/2) 



-s 



T{u)T(k + (l/2)di)«i 

^ (fcp + 1/2), / 2fa + 2£+dA a / /l \ 

X ^ (^(l72KM27)! r [—2^ ) ^ ^"(^M (,2 s ; (32) 

1 ^ (fc +£)!(l/2) fco+ ,c fco+ ,( ai ) ^ 2 £ /l \ 

^f(^^ (2^! ( x) ^'-^/^{r) 



oo 



1 ^ m! (V 2 )mC m (ai) , 2 ^ /" 1 \ 

» 2^ {2(m-fc )}! ( ^^'-"V"!^- 



«ir(i/) 

Next consider G™(0,y) for 7n < fco. The process |^Z(x, y) has spectral density 
r 2m (|x| 2Ql +|a;| 2 )- ,y , and hence the process -§£nZ(0, y) = Zf(0,y) considered just as a 
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dl (| T |a«i + | w |2)« 

2 7I (di-l)/2 roo |-7t w 2m+di-l cos 2m ^ 

= r((d! - 1)/2) i /„ 

7t (rfi- 1 )/ 2 r(m + l/2)B((2m + d 1 )/2ai,z/-(2m + di)/2Q!i) 



air(m+ (l/2)<ii)|u;| 2 ^-( 2 ™+ d i)/ Q i 

by switching r to hyperspherical coordinates and using [8], 3.241.4. The process Z™(Q,y) 
with spectral density (33) is an IRF-(fco — m) (it may be an IRF of lower order as well) 
and its corresponding GC-(fco — m) can be taken as [6], Chapter II, Section 3.3, equations 
(2) and (11), 



G™(0 y) = pm(Q!x)m!r(m + l/2) 



l9'-,n/ ai (l}\y\j 



n 1 / 2 a 1 T(iy) 

(34) 

c m ( ai )m!(l/2) m (1 . 

= T^'-W^M 

To recover Gi(xi,y) and hence G(x, y), repeatedly integrate (32) and use (34) to set 
the boundary conditions. Specifically, Gi(xi,y) must be of the form 

Gi(n,y)=Wf<(y) + (-i) h f 1 I" 1 -- p k ° 1 G k 1 °(z 2ko ,y)dz 2ko ---dz 1 

e=Q J<3 J<3 JO 

for some suitable functions i*b, . . . ,-Ffc _i. Substituting the series for G\° in (32) into 
the preceding expression and integrating termwise, which is easily justified for s > by 
dominated convergence, yields 

fco-l 

Gi(si,y) = J>fF/(y) 

(35) 

1 c m (a 1 )m!(l/2) m 2 2 "' f /l \ 2 1™ (1 A 



Elementary calculations demonstrate 
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for all m € No- For < m < fco, differentiating (35) 2m times, setting y = and using 
(34) and (36) gives 



_ (-l/4) m c m ( ai ) 



aiT(u) 

Substituting (36) and (37) into (35) yields 



(37) 



Gi(xi,y) = Yl {^Q^ 1 )} c "^ a ^s'-m/ ai Q|y 

and (10) follows from Gi(xi,y) = G(|aii|, |y|). 

To obtain an explicit expression for G(r, 0), go back to (25) and change the order of 
integration. Integrating by parts and using [8], 3.761.4, 

XlCOS<^l f°° 2An-2r>,fl' • / , \ i 
r / u zl ZaiV sm (a; 1 MCOS0i)dM 



2ai0' - 2fc Jo 
= -r(2fc - 2a x 6') sin{7t(fc - ax&)}\xi cos fa \ 2a ^'- 2k '> . 

When fc ^ ot\0' , substituting this result and M.e(Q) = 2 e ~ 1 T(8) into (25) and integrating 
over 0i yields 

lfco _ 27t^+^- 1 )/ 2 r(0)r(a 1 0' + l/2)r(2fc - 2a i / ) 



Gf (x ls 0) 



r(i/)r(aifl) 

x cos(7ta 1 6'')(-l)' ;o |xi| 2Qie '- 2fco . 



Now, r(2fco — 2ai0') = (2ai0' — 2fco + l)2fc o r(— 2a\6'), so applying the duplication formula 
for r to T(—2a.i6') and then the reflection formula to T{—ot\Q' + ^) yields 

Gj°(zi,0) 

_ (-l) fc °7r( rf 1 +^)/ 2 r(g)r(-a 1 0O(2 Ql ^ - 2fc + l) 2fco 2ai0 ,_ 2feo 

r(i/)r( ai 6i)2 Q i e ' 1 11 

Integrating this expression 2fco times and using the boundary condition G™(0, 0) = for 
m < fco to make G™(0,y) continuous at y = gives (11). 

To show that (10) holds when fco = ol\B' , write G„ to make the dependence of 
G on v explicit (but still suppressing the dependence on ai, d\ and d 2 ). Define 
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v o = (ko + ^d\)/ot\ + \d^ and view ct\ as fixed. For any given ALC-fco, the last line 
of (9) is continuous as v \, i>o, so the first line is as well. Thus, 

n n 

V \ i \jGv (rej,st j )= lim V A^A i (5 I ,(r«,s^)- ( 38 ) 

' — ' viva * — ' 

fj=i e,j=i 

Let Mo be the set of nonnegative integers m for which h m = (ko —m)/a\ e No. This set 
is finite and includes fco as its largest element. Writing v — i/q + e, define 

Because y? —1 A^Z(x^) is an ALC-fco, for a, & £ No and a + b < fco, subtracting any linear 
combination of terms like rf"s|* from G v (rgj,S£j) on the right-hand side of (38) does 
not change the result. Now ot\ > 1 implies m + h m < fco, so P e (r£j, sij) is of this required 
form and 

n n 

Thus, to prove (10), it suffices to show that for all r,s nonnegative, 

G, (r, s) = lim{G„(r, s) - P e (r, «)}. (39) 

For s > 0, for all £ sufficiently small, 

2n rn 



^.•>-'.<'.->-£^{-(H } 



( r s 



2/i„ 



2f 



2 S ' 



Dominated convergence justifies taking the limit v i/q (equivalently, as e J, 0) inside 
this infinite sum. For m ^ Mo, the limit is trivial, so consider m G Mo- By the reflection 
formula for T, for n e No, 

7r(-l) ,l+1 

r(-n-e) = — , , r. (40) 

V ' sin(7te)r(n+l + e) V ; 

Using this result and straightforward calculus yields 

^ \ 2h m ( / ]\ 2e "1 / ^ 



as v I Vq , establishing (39) when s > 0. 
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It remains to establish (39) when s = with G Va (r, 0) given by (12). For all e sufficiently 
small, applying (40) to T(— fco — otie) and using T(l — e) = — eT(— e), for all e sufficiently 
small, 



G„(r,0)-P s (r,0) 

( _ 1) fco+l 7T (d 1 +d 2 +2)/2 r ^ :) 

~ sin(7tai£r)r(i/)r(aifl 1/ )r(fc + 1 + a x e) \ 2 



2k +2a 1 e 



_ c ko {ai)T{-e) 

n (d 1 +d 2 )/2 



2 \ ko 



2-sk 



(41) 



nT{0 v ) 



—r 



sm(nais)T(ai9 u )r(k + 1 + cue) l{\2 
r((d 1 + 2ft )/2a 1 )r(l-e) 



2aie 



-i+i 



aier(fc + (l/2)d 1 )A:o! 



As e|0, 



71 



1 



sin(7taie) [ \ 2 



- 1 } -»21og| -r 



(42) 



By the definition of the digamma function tp, T(x + e) — r(x){l + etp(x) + 0(e 2 )} as 
e — s- as long as — x £ Nq. Then 



7rr(0„) 



r((di+2*o)/(2ai))r(l-e) 



sin(7taie)r(ai6' ly )r(fco + 1 + ens) aieT(k + (l/2)d 1 )k \ 
_ r((di + 2fco)/(2ai)) 



l-eV(l) 



+ 0(e) (43) 



fc !r(fc + (l/2)di) 

x [ l+eV;((di+2fco)/(2ai)) 

[aie{l + aie^Cfco + (l/2)di)}{l + aie-0(fc o + 1)} cue 

= r((di+2fc )/(2ai)) 
fc !r(fc + (l/2)dr) 

x + ^V'(l) -1>(ko + \d^j 4>{k, + 1) j + 0(e). 

Thus, when fco = ai^' and G Vo {r, 0) is defined as in (12), (41)-(43) imply (39) holds for 
s = 0. 
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A. 2. Properties of ii-functions 

This Appendix provides some background material on i?-functions and is taken 
from [14], Section 1.1. Suppose m, n,p and q are integers satisfying 0<m<g,0<n<p, 
ai, . . . , a p , b\, . . . b q are complex numbers and ai,...,a p ,/3i,...,/3 g are positive reals. Then 
the ff-function Hi 71 '™ is defined by, for complex z, 



rrm^n 



(oi,ai), . . . , (a p ,a p ) 
ft),..., 

2m Jc iHU+i r(a, + M IlLm+i r ( X " &i + M 



(44) 



; ds, 



where [14], p. 2, gives the form of the contour £, and an empty product is defined to be 
1. For this integral to be well defined, none of the poles of the gamma functions in the 
two products in the numerator of (44) may coincide, or 

<*.,(&,+ ^ ft K-fc-1) (45) 

for 1 < i < n, 1 < j < m and all k,£ £ Nq. The validity of series expansions of H™^ 1 
generally depends on the signs of the following two quantities: 

n p m q 

a *=E a '*- E «*+Eft- E ft ( 46 ) 

i—l i— n+1 i=l j=m-|-l 

and 

A =Eft-E a - ( 47 ) 

j=l t=l 

where an empty sum is defined to be 0. 
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